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AN ALGEBRAIC MULTIGRID METHOD FOR 
QUADRATIC FINITE ELEMENT EQUATIONS OF 
ELLIPTIC AND SADDLE POINT SYSTEMS IN 3D 

HUIDONG YANG 


Abstract. In this work, we propose a robust and easily imple¬ 
mented algebraic multigrid method as a stand-alone solver or a pre¬ 
conditioner in Krylov subspace methods for solving either symmet¬ 
ric and positive definite or saddle point linear systems of equations 
arising from the finite element discretization of the vector Lapla- 
cian problem, linear elasticity problem in pure displacement and 
mixed displacement-pressure form, and Stokes problem in mixed 
velocity-pressure form in 3D, respectively. We use hierarchical qua¬ 
dratic basis functions to construct the finite element spaces. A 
new heuristic algebraic coarsening strategy is introduced for con¬ 
struction of the hierarchical coarse system matrices. We focus on 
numerical study of the mesh-independence robustness of the alge¬ 
braic multigrid and the algebraic multigrid preconditioned Krylov 
subspace methods. 


1. Introduction 

Compared to the geometrical multigrid (GMG) method (see, .e.g., 
i), the algebraic multigrid (AMG) method (see, e.g., [HI [TS]) is a 
purely matrix-based approach, that does not rely on any underlying 
mesh hierarchy; see, e.g., |H] for the development from GMG to AMG 
methods. Goncerning comparison of different types of AMG methods 
we refer to, e.g., 122 ] for a review and related references. In contrast 
to coarsening based on the strongly connected matrix entries in the 
classical AMG method, an AMG method (among others) with special 
coarsening and interpolation strategies was introduced in HH, that 
is based on graph connectivity of the matrix only and leads to fast 
construction of matrices on coarse levels. An AMG method that is 
based on the matrix graph information only was also studied early in 
[3]. The further development of such an AMG method [H] in different 
applications have been reported in, e.g., Ha ESI ESI [la [H [29]. In 

Key words and phrases, algebraic multigrid method, algebraic multigrid precon¬ 
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this work, we focus on the development of such an AMG method for 
both elliptic and saddle point systems of equations arising from the 
quadratic hnite element discretization for the three dimensional (3D) 
vector Laplacian problem, linear elasticity problem in pure displace¬ 
ment and mixed displacement-pressure form, and the Stokes problem 
in mixed velocity-pressure form. This requires new coarsening strate¬ 
gies to construct the hierarchy of matrices on coarse levels for both 
the elliptic and saddle point systems, that are to be developed in this 
work. We notice, that different AMG methods towards higher-order 
hnite element equations for second order elliptic problems were also 
studied using different approaches in, e.g., m [16]. The main focus 
of this work is the numerical study of the robustness and efficiency of 
the designed AMG method as a stand-alone solver or a preconditioner 
in Krylov subspace methods for solving the elliptic and saddle point 
systems. 

The remainder of this paper is organized in the following way. In 
Section we describe the model problems, their hnite element dis¬ 
cretizations and the arising linear systems of equations. The alge¬ 
braic multigrid method using a new heuristic coarsening strategy is 
prescribed in Section In Section we present numerical results of 
the AMG method applied to discrete model problems. Finally, some 
conclusions are drawn in Section |5l 

2. Preliminaries 

2.1. The model problems. Let G C be a simply connected and 
bounded domain with two boundaries Tjv and T ^ such that T^) U f tv = 
dil and TTvnT d = 0- We consider the 3D vector Laplacian problem, the 
linear elasticity problem in pure displacement and mixed displacement- 
pressure forms, and the Stokes problem in mixed velocity-pressure 
form, that are formulated in the following: 

For the vector Laplacian problem: Find the potential u : D i—>■ 
such that 

(1) — Au = 0 in D 

with the boundary conditions u = gD ouTd and |^ = on Ftv, where 
n denotes the outward normal vector on Ftv- 

For the linear elasticity problem in pure displacement form: Find 
the displacement u : D i—)■ such that 

(2) — V ■ a{u) = 0 in D 

with the boundary conditions u = gr, on Vr) and a{u)n = ^fTv on Ftv- 
In particular, we use the linear Saint Venant-Krichoff elasticity model. 
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The Cauchy stress tensor and the inhnitesimal strain tensor are defined 
by a{u) = 2fj,e{u) + Adiv(M)/ and e{u) = (Vu + Vm^)/ 2, respectively, 
with Lame constants A and fi. 

For the linear elasticity problem in mixed displacement-pressure form: 
Find the displacement u : i—)■ and pressure p : 12 i—)■ M such that 


(3) 


-V ■ {2fj,e{u)) -|- Vp = 0 in 12 


-V ■ M — —p = 0 in 12 
A 


with the boundary conditions u = qd on Td and {2^e{u) — pl)n = 
on Ftv- It is easy to see that, in this classical mixed displacement- 
pressure form, the displacement and pressure are associated by the 
relation p = —AV ■ u; see, e.g., |1]. 

For the Stokes problem in mixed velocity-pressure form: Find the 
velocity u : 12 i—>■ and pressure p : 12 i—)■ M such that 


(4) 


— V ■ {2pe{u)) -I- Vp = 0 in 12 
—V ■ n = 0 in 12 


with the boundary conditions u = qd on To and {2pe{u) — pl)n = 
on Ftv, where p denotes the dynamic viscosity. 


2.2. The variational formulations. We search for weak solutions 
of the above four model problems ([^-Q in proper spaces. For this, 
let iL^(12) and Q = L^(12) denote the standard Sobolev and Lebesgue 
spaces on 12; see m- With V = Ff^(12)^, we define the spaces = {u G 
V : mIto = go} ioT the potential and displacement (velocity) functions. 
We also define the homogenized space Vq = {u & V ■. u\yj^ = Q}. In 
addition, we assume the given data go G d) , where d) 

denotes the trace space, i.e., = {n|ro • ^ iL^(12)^}. We 

also assume the given data gN G L^(Fat)^. By standard techniques, the 
following variational formulations are obtained. 

The variational formulation for the vector Laplacian problem ([^ and 
the linear elasticity problem ([^ in pure displacement form reads (after 
homogenization): Find u E Vq such that 

(5) a{u,v) = {F,v) 

for all V G Vq, with the bilinear form a{u,v) := f^Vu : Vvdx for the 
vector Laplacian problem and a{u,v) := J^[2pe{u) : e(v) + AV ■ uV ■ 
v\dx for the linear elasticity problem, respectively, and the linear form 
(F, n) := J-p^gN ■ vdx — a{gD,v), accordingly. 

The variational formulation for the linear elasticity problem (|^ in 
mixed displacement-pressure form and the Stokes problem (|^ in mixed 
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velocity-pressure form reads (after homogenization): Find u G Vq and 
p ^ Q such that 

a{u,v)+ b{v,p) = {F,v), 
b{u,q) - c{p,q) = {G,q) 

for all V E Vo and q E Q, where the bilinear and linear forms are 
given by a{u, v) = 2p, e(u) : e(v)dx, b(v, q) = — Jq ■ vdx, {F, v) = 
gjv-vdx — algBjV), {G,q) = 0 , respectively, and c{p,q) = j f^pqdx 
and c{p, q) = 0 for the linear elasticity and Stokes problem, respectively. 

2.3. The finite element discretization. The spatial discretization 
is done by the Galerkin hnite element method with a hierarchical qua¬ 
dratic polynomial basis functions. Let Fh be the admissible subdivision 
of the domain hi into tetrahedra. The four linear basis functions on each 
tetrahedron T E Th are nothing but standard Pi hat functions in 3D, 
i.e., (pi := Aj, i = 1,...,4, where A* are the barycentric coordinates of 
T. The six qnadratic basis fnnctions are then dehned as 

05 = 4 A 1 A 2 , 06 = 4 A 2 A 3 , 07 = 4 A 3 A 4 , 

08 = 4 A 1 A 3 , 09 = 4 A 1 A 4 , 010 = 4 A 2 A 4 , 

that constrnct hierarchical qnadratic polynomial basis fnnctions. 

Let Vl ■.= {v E Go(D) : vt E <Fr,VT G Th\ be the snbspace of 
continnous piecewise linear hat fnnctions with zero traces on Td and 
Vq := {n G Go(D) : vt E E Th} the snbspace of continnous 

piecewise quadratic functions with zero traces on Td, where $ 7 ’ = 
span{ 0 i=i ,...,4 : T E %] and Tr = span{ 0 i= 5 ,„.,io : T E %}■ Let 
Ql ■= e G(f2) : Vt E <Fr,VT G Th} be the subspace of continuous 
piecewise linear hat functions. It can be shown that the global degrees 
of freedom (DOF) of Vl, Vq, Ql are the nnmber of vertices (/) and 
edges (m) (exclnding the vertices and edges on F^), and the nnmber 
of all vertices (n), respectively. The global basis fnnctions ipi E Vl,Ql 
and tpi E Vq can be constructed from the local ones. The function 
space for one component of the potential or displacement is dehned as 
Vh = Vl ® Vq E Vq, a. linear snbspace complemented by a qnadratic 
snbspace. The fnnction space for the pressnre is dehned as Qh = Ql F 
Q. 

Using Galerkin’s principle the discrete elliptic variational formulation 
for the vector Laplacian and linear elasticity problem in pure displace¬ 
ment form read: Find Uh E 14 snch that 

(7) a { uh , Vh ) = { F , Vh ) 

for all Vh E Vh- 
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The discrete mixed variational formulation for the elasticity problem 
in mixed displacement-pressure form and the Stokes problem in mixed 
velocity-pressure form reads: Find {uh,Ph) ^VhXQh such that 

a{uh,Vh) + b{vh,Ph) = {F,Vh), 

Qh) - c{ph, Qh) = {G, Qh) 

for all Vh e Vh and qn^Qh- 

The hnite element solutions Uh G 14 and ph € Qh are expressed by 
the ansatz: 

I m n 

Uh = u{ + ul = ^u\pi + ^ul'4)i, Ph = '^p\^i, 
i=l i=l i=l 

respectively, where u\,ul G and p\ G M. It is easy to see the hnite 
element solution Uh is the sum of the linear and quadratic part, 
and respectively. For the pressure ph, we have a linear approxima¬ 
tion. The mixed hnite element for the elasticity and Stokes problem 
is classical Taylor-Hood element, that fulhlls the inf — sup stability re¬ 
quirement; see, e.g., [6]. 


2.4. SPD and saddle point linear systems of eqnations. Using 
the hnite element discretization (including homogenization), we obtain 
the following symmetric and positive dehnite (SPD) system of equa¬ 
tions for the elliptic problem: 


(9) 


Ku 



1 

IS 


1- 

'4T 

1_ 

Kgi 
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LhJ 


/, 


where Ku = Kgi = a{ipj,'il!i), Kgg = Ui = («•), 

Mg = h = {{F,Pi)) and = ((F,^^)). 

For the elasticity problem in mixed displacement-pressure form and 
the Stokes problem in mixed velocity-pressure form, we obtain the fol¬ 
lowing symmetric indehnite system of equations: 

( 10 ) 


r A 

B^ 1 


B 

-C 






"k 



Ku 


Bfi 

Kgi 

K,, 

U 

Bu 

Big 

-Cii 


' M/ ' 


1 


= 

K 

aj 


L h J 


/ 

9 


where Ku = Kgi = Kgg = Bu = 

6((pi,(pj), Big = b{ipi,ijj), Cii = c{ipi,iph), m = {u{), Ug = (uj), 

h = Li = Lg = and = {{G,p^)). It is 

obvious that for the Stokes problem, (7 = 0. 
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In the following, we focus on how to solve the above two systems of 
equations ([^ and (10) using an algebraic multigrid method. 


3. An algebraic multigrid method 

3.1. The basic algebraic multigrid iteration. The basic AMG it¬ 
eration applied to a general linear system of equations Kx = b is given 
in Algorithmic with rupre and rripost being the number of pre- and post¬ 
smoothing steps (steps 1-3 and 14-16, respectively). By choosing u = 1 
and V = 2, the iterations in Algorithmic are called V- and W-cycle, re¬ 
spectively. As a convention, we use / = 0,..., L to indicate the algebraic 
multigrid levels from the hnest level / = 0 to the coarsest level I = L. 
On the coarsest level L, the system is solved by any direct solver (step 
6). The coarse grid correction step is indicated in steps 4-13. The full 
AMG iterations are realized by repeated application of this algorithm. 
The iteration in this algorithm is also combined with the Krylov sub¬ 
space methods, that usually leads to accelerated convergence of V-cycle 
or W-cycle preconditioned methods [22] . 


Algorithm 1 Basic AMG iteration: AMG(A';, 6;) 
I: for k = 1 to rUpre do 
2 : = Si{x^i ,hi) 

3: end for 

4: hi+, = R\^\hi-KiXi), 

5: if 1-|-1=L then 
6: Solve KlXl = 

7: else 

8: Xz+i = 0, 

9: for fc = 1,..., z/ do 

10: xi+i =AMG{Ki+i,xi+i,bi+i), 

11: end for 

12: end if 

13: Xi=Xi + P/^iTz+i, 

14: for = 1 to mpost do 
15: xf^^ = Si{x’l,bi), 

16: end for 

17: return Xi. 


3.2. A new heuristic coarsening strategy. 
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3.2.1. Case I : The SPD system. A robust coarsening strategy is an 
important feature of the AMG method, that is used to construct the 
system matrices on coarse levels I = For the second order 

elliptic equations discretized by low order hnite element or boundary 
element methods, some well known graph-based black-box or grey- 
box type AMG methods have been introduced and applied, see, e.g., 
[miElGSlE]. The general strategy is to split the nodes into the sets 
of coarse and hne nodes, based on the graph connectivity of the system 
matrix or the constructed auxiliary matrix (’’virtual” hnite element 
mesh [T7]). 

When applying such a technique to the system matrix in ([^, we 
obtain a very dense graph connectivity constructed from the stiffness 
matrix A, that contains connectivities for the linear DOF, the quadratic 
DOF and the coupling between them. This will lead to a mixture of 
different oder of DOF and may cause additional difficulty to construct 
the interpolation operators. In fact, from our numerical studies, we 
observe the loss of optimality of the AMG method when such a dense 
graph connectivity is adopted for coarse system matrix construction. 
Therefore, we construct the graph connectivities only for the linear and 
quadratic DOF, i.e., for Ku and Kgg, respectively. We simply neglect 
the coupling connectivity in Kqi. By this means, we avoid the mixture 
of different order of DOF on the coarse level. In addition, we are able 
to construct and control the interpolation operators for the linear and 
quadratic part, respectively. A simple comparison of the classical and 
new graph connectivities is illustrated in Fig. 



Figure 1. Graph connectivity constructed using the 
new (left) and the classical (right) strategies: linear DOF 
(solid line) , quadratic DOF (dashed hne) , coupling 
(dashed dot lines). 
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From the left plot in Fig. [TJ we show the two graph connectivities 
indicated by solid and dashed lines for the linear and quadratic DOF, 
respectively. For a comparison, on the right plot, the graph connectiv¬ 
ities constructed by the new strategy is reconstructed by the classical 
strategy. For simplicity, we only reconstruct the part indicated by the 
lines with blue color. It is easy to obverse that, the classical one leads to 
much denser graph connectivities than the new one due to the coupling. 

Based on these two graph connectivities, the prolongation matrix 
Pij^i from the coarse level / -|- 1 to the next hner level I is constructed 
in form of 

(11) A+i = 

where the prolongation matrices, —)■ and : 

(]^3)m;+i are dehned for the linear and quadratic DOF, re¬ 

spectively, where ni and mi denote the number of linear and quadratic 
DOF on level /, respectively. The restriction matrix from the hner level 
I to the next coarser level / -|- 1 is constructed as . The system 

matrix Ai+i on the level / -|-1 is constructed by the Galerkin projection 
method: 

We mention that the two graph connectivities for the linear and qua¬ 
dratic DOF are naturally different as illustrated in Fig. that may 
require different coarse and hne nodes selection algorithms. However, 
for simplicity, we apply the coarse and hne nodes selection and prolon¬ 
gation operator matrix construction algorithms developed in m for 
both the linear and quadratic DOF, that show the robustness from the 
numerical studies. 

3.2.2. Case II: The saddle point system. A robust coarsening strategy 
for saddle point problems is, in general, more involved than that for 
elliptic problems mainly due to the inf — sup instability issues possibly 
caused by standard Galerkin projection method. For saddle point prob¬ 
lems arising from the low order hnite element discretized huid problem, 
the stability issue has been studied in, e.g., [2211251 [U] . [211ES], a so- 

called 2-shift coarsening strategy was introduced for the discrete huid 
problem using the modihed Taylor-Hood element (PiisoP 2 — Pi), that 
mimics the hierarchy of matrices in the geometrical multigrid method. 
To guarantee the inf — sup stability for the coarse system is still a re¬ 
search topic. 

We extend the new coarsening strategy described above for the el¬ 
liptic problem, to the saddle point problem, based on the new graph 


r 


^l+l 
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connectivity construction. As illustrated in Fig. we show the graph 
connectivities constructed for the velocity (displacement) and pressure. 
For the velocity, we follow the same strategy as for the SPD system; see 
the left plot. For the pressure, we have conventional graph connectivity 
for the linear hnite element matrix; see the right plot. 



Figure 2. Graph connectivities constructed using the 
new strategies for velocity (displacement) (left) and pres¬ 
sure (right). 


Based on the graph connectivities, the prolongation matrix from 
the coarse level / -|- 1 to the next hner level I is constructed in form of 

(12) Gi = 




Ji 


/+! 


L J 

where the prolongation matrices, : (M^)"''+i —)■ (M^)""* and : 

are dehned for the linear and quadratic velocity 
DOF, respectively, where rii and mi denote the number of linear and 
quadratic velocity DOF on level I, respectively, —)■ for 

the pressure DOF, where ki the number of pressure DOF on level 1. 
The restriction matrix from the hner level I to the next coarser level 
/ 1 is constructed as The system matrix on the level 

/ 1 is constructed by the Galerkin projection method: 


A-,+1 = (A+i)’’AT’/+i. 


We admit that the inf — sup stability of the coarse system is still open 
by this construction. However, from the numerical studies, we observe 
quite satisfactory results using this new coarsening strategy. Neverthe¬ 
less, we are at least able to obtain efficient multigrid preconditioners 
by using pure Galerkin projection; see comments in, e.g., [211 [2Z] and 
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numerical experiments for the Stokes problem using the low order hnite 
element discretization in, e.g., [10]. 


3.3. The smoothing procedure. To complete the algebraic multi¬ 
grid algorithm, a smoothing procedure is needed. As conventional 
choices, we employ the damped block Jacobi and block Gauss-Seidel 
smoothers for the SPD system ([^, that are widely used in the multi¬ 
grid methods. For the saddle point system (10), we have considered 
the following smoothers, that were originally designed and analyzed in 
the GMG method. 


3.3.1. The multiplicative Vanka smoother. The multiplicative Vanka 
smoother was introduced in [23| for the fluid problem. We have recently 
developed an AMG method with this smoother for solving the nonlinear 
and nearly incompressible hyperelastic models in fluid-structure inter¬ 
action simulation [IT]. To adapt this smoother for the Taylor-Hood 
element, we hrst construct the patches Vi, i = l,...,n. Each patch 
contains one pressure DOF, and the connected linear and quadratic 
velocity DOF indicated by the connectivity of matrices Bu and Big, 
respectively. A typical patch Vi is illustrated in Fig. 



pressure DOF 

O linear velocity DOF 
(') quadratic velocity DOF 


Figure 3. A typical local patch Vi for the Taylor-Hood 
element contains one pressure DOF, and connected linear 
and quadratic velocity DOF. 
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The local (correction) problem on Vi is extracted by a canonical 
projection of the global one (10) to local one on Vi. 


(13) 


fc +1 

'i 

k-\-l 


u- 

k 

Pi 


+ OJ 




Bf 

-Ci 


-1 r 


' p,i 


with k representing the smoothing step and oj being a damping pa¬ 
rameter. Here [(r^ denotes the residual updated in a mul¬ 

tiplicative manner. As we observe from the numerical studies, this 
smoother shows the efficiency and robustness if it is used in the AMG 
preconditioner but not in the stand-alone AMG solver. 


3.3.2. The Braess-Sarazin-type smoother. This smoother was introduced 
in |5] and approximated in [30], that has been applied to the fluid prob¬ 
lem [25l [26] , the nearly incompressible elasticity problem [28] and the 
fluid-structure interaction problem [22] • One smoothing step corre¬ 
sponds to a preconditioned Richardson method: 


(14) 


fc+i 

fc+i 


u 

k 

P 


+ K 


-1 


/ — — B^p^ 

g — Bu^ + Cp’^ 


with preconditioner 


(15) 


k 


A B^ 

B bA-^B^-S 


As in 129, we use A = 2D, where D represents the diagonal of A. We 
use an AMG preconditioner S for the approximated Schur complement 

c + bA-^b. 


3.3.3. The segregated Gauss-Seidel smoother. This smoother was very 
recently introduced in [7] as a segregated Gauss-Seidel smoother based 
on a Uzawa-type iteration. Such a Uzawa method (see, e.g., 0) can 
be reinterpreted as a preconditioned Richardson method: 


(16) 


fc+l 

fc +1 


u 

k 

P 


+ K 


-1 


/ — Ak — B'^p^ 
g — Bk + Cp'^ 


with the preconditioner 


(17) 


k 


A 

B -00-^1 ’ 


where a; is a properly chosen parameter, A is a (e.g., AMG) precon¬ 
ditioner for A. However, to get a multigrid smoother, this is relaxed 
by choosing some proper smoother for A instead of a precondi¬ 
tioner A; see [7] . The theoretical analysis requirement for has been 
specihed therein. In our setting, we have chosen Ma as a damped 
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block Jacobi smoother with a damping parameter 0.5. Compared to 
the Braess-Sarazin-type smoother, this smoother avoids the explicite 
construction of the approximate Schur complement. 

4. Numerical results 

4.1. Meshes, boundary conditions and coarsening for the vec¬ 
tor Laplacian and linear elasticity problems. We consider a unit 
cube (0,1)^ as the computational domain for the vector Laplacian and 
linear elasticity problems. The domain is subdivided into tetrahedra 
with four levels of mesh rehnement Li — L 4 ^. The number of tetrahe¬ 
dron (#Tet), nodes (#Nodes) and midside nodes (^^Midside nodes), 
and the total number of DOF for elliptic ( (elliptic) ) and sad¬ 

dle point ((#DOF (saddle point)) systems are shown in Table We 


Level 

Li 

L 2 

L 3 

T 4 

#Tet 

64 

512 

4096 

32768 

#Nodes 

125 

729 

4913 

35937 

^(/^Midside nodes 

604 

4184 

31024 

238688 

#DOF (elliptic) 

2187 

14739 

107811 

823875 

#DOF (saddle point) 

2312 

15468 

112724 

859812 

Table 1. Number of 

tetra 

ledron 

(#Tet) 

, nodes 


(#Nodes), Midside nodes (^Midside nodes), and total 
number of DOF for elliptic ( ^DOF (elliptic) ) and sad¬ 
dle point (#DOF (saddle point)) systems on four levels 
Li — L4. 


£x the bottom of the domain, i.e., u = [0,0,0]^ at 2 ; = 0, prescribe a 
Dirichlet data on the top, i.e., u = [0, 0,1]^ at 2 ; = 1, and use zero Neu¬ 
mann condition on the rest of the boundaries. For the linear elasticity 
problem, we set p = 1.15e -|- 06 and A = 1.73e -|- 06. In Fig. we plot 
the value of ||m||r3 of the numerical solution u (indicated by the color) 
for the vector Laplacian (left) and the linear elasticity problem in pure 
displacement form (right), respectively. For visualization purpose, we 
plot the vector helds of potential and deformation, that is scaled by a 
factor of 0.1. 

As a comparison, on each algebraically coarsening level, we show the 
number of linear and quadratic DOF (^ Linear DOF and Quadratic 
DOF, respectively) in the new coarsening strategy, and the classical 
coarsening strategy (# Non-separating DOF); see in Tablethe num¬ 
ber of DOF on each coarsening level for the vector Laplacian and linear 
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Figure 4. Numerical results of the vector Laplacian 
(left) and linear elasticity (right) problem. 

elasticity problem on the level L 4 . It is obvious to see these two strate¬ 
gies lead to different graph connectivity in the coarsening procedure. 


Goarsening Levels 

0 

1 

2 

3 

4 

^(/^Linear DOF 

107811 

14739 

2187 

375 

81 

^Quadratic DOF 

716064 

104544 

6366 

369 

39 

^(/^Non-separating DOF 

823875 

88467 

2187 

375 

24 


Table 2. Number of linear and quadratic DOF in the 
new coarsening strategy, and the total DOF in the non¬ 
separating coarsening strategy at level L 4 , for the vector 
Laplacian and linear elasticity problems. 


4.2. Numerical performance for the vector Laplacian prob¬ 
lem. Before showing the AMG performance with the new coarsen¬ 
ing strategy, we demonstrate the performance with a black-box type 
AMG m in Table for the vector Laplacian problem, where non¬ 
separating coarsening strategy is used. For both AMG and AMG 
preconditioned GG methods, we use the relative residual error ||/ — 
^'^kWh/Wf ~ ~ Ih® / 2 —norm as stopping criteria, 

where k denotes the number of AMG iterations. In each iteration of 
the AMG solver, we use 2 W-cycles and 1 pre- and post-smoothing 
step. As observed, the AMG and the AMG preconditioned GG are not 
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robust with respect to the mesh rehnement, i.e., the iteration number 
increases with mesh rehnement. 


Levels 

Li 

L 2 

L 3 

F 4 

#It AMG 

90 

158 

> 200 

> 200 

#It PGG AMG 

22 

25 

36 

64 


Table 3. Performance of black-box type AMG solver 
(^^It AMG) and AMG preconditioned conjugate gradient 
solver (^^^It PGG_AMG) for the vector Laplacian prob¬ 
lem. 


Now we show the performance of the AMG and AMG precondi¬ 
tioned GG solvers using the new coarsening strategy. Note, that in 
the following numerical tests, we stop the iterations when the relative 
residual error in the ^ 2 —norm is reduced by a factor 10^^. We consider 
the Jacobi smoother with damping parameter 0.5 and 1 or 2 pre- and 
post-smoothing steps (JA-1-1-0.5 or JA-2-2-0.5), and the Gauss-Seidel 
smoother with 1 or 2 pre- and post-smoothing steps (GS-1-1 or GS- 2 - 
2). In Table we show the performance of the AMG solver for the 
vector Laplacian problem using V-cycle with different smoothers. In 
Table we show the performance using W-cycle. In Table and we 
show the performance of the AMG preconditioned GG using V- and 
W-cycles with different smoothers, respectively. 

As observed, the iterations for each solver are independent of mesh 
rehnement levels. The AMG solver using the Gauss-Seidel smoother 
shows better performance than the damped Jacobi smoother. By using 
the GG acceleration, we observe similar performance with two differ¬ 
ent smoothers. In addition, we observe, that the V- and W-cycles 
demonstrate almost the same performance. We also observe that the 
computational cost is proportional to the number of DOF. 


Table 


Level 

Li 

L 2 

L 3 


#It ( JA-1-1-0.5 ) 

129 

125 

124 

128 

#It ( JA-2-2-0.5 ) 

65 

65 

65 

66 

#It ( GS-1-1 ) 

43 

46 

47 

47 

#It ( GS-2-2 ) 

23 

24 

24 

24 

4. Performance 0 

■ the 

AMG so 

ver ; 


or the 


vector Laplacian problem using V-cycle with different 
smoothers. 
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Table 


Level 

Li 

L 2 

L 3 

U 

#It ( JA-1-1-0.5 ) 

128 

124 

121 

123 

#It ( JA-2-2-0.5 ) 

65 

64 

63 

64 

#It ( GS-1-1 ) 

43 

46 

47 

47 

#It ( GS-2-2 ) 

23 

24 

24 

24 

5. Performance o 

■ the 

AMG so 

ver ; 


'or the 


vector Laplacian problem using W-cycle with different 
smoothers. 


Level 

Li 

L 2 

L 3 

L, 

#It ( JA-1-1-0.5 ) 

30 

30 

30 

30 

#It ( JA-2-2-0.5 ) 

21 

22 

22 

22 

#It ( GS-1-1 ) 

26 

29 

29 

30 

#It ( GS-2-2 ) 

17 

19 

19 

19 


Table 6. Performance of the AMG preconditioned CG 
solver for the vector Laplacian problem using V-cycle 
with different smoothers. 


Level 

L, 

L 2 

L 3 

L, 

#It ( JA-1-1-0.5 ) 

30 

30 

30 

29 

#It ( JA-2-2-0.5 ) 

21 

21 

21 

21 

#It ( GS-1-1 ) 

26 

29 

29 

29 

#It ( GS-2-2 ) 

17 

18 

18 

18 


Table 7. Performance of the AMG preconditioned GG 
solver for the vector Laplacian problem using W-cycle 
with different smoothers. 


4.3. Numerical performance for the linear elasticity problem 
in pure displacement form. We perform the same test for the lin¬ 
ear elasticity problem. In Table we show the performance of the 
AMG solver for the linear elasticity problem using V-cycle with differ¬ 
ent smoothers. In Table we show the performance using W-cycle. In 
Table IT and[T^ we show the performance of the AMG preconditioned 
GG using V- and W-cycles with different smoothers, respectively. 

As observed, the damped Jacobi smoother does not work for this 
test problem. The AMG solver using the Gauss-Seidel smoother shows 
good performance. By using the GG acceleration, we observe improved 
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performance. We observe, that the V- and W-cycles demonstrate al¬ 
most the same performance. 


Level 

Li 

L 2 

L 3 


#It ( JA-1-1-0.5 ) 

— 

— 

— 

— 

#It ( JA-2-2-0.5 ) 

— 

— 

— 

— 

#It ( GS-1-1 ) 

80 

78 

76 

75 

#It ( GS-2-2 ) 

44 

40 

39 

39 


Table 8. Performance of the AMG solver for the linear 
elasticity problem in pure displacement form using V- 
cycle with different smoothers.. 


Level 

L, 

L 2 

L 3 

L, 

#It ( JA-1-1-0.5 ) 

— 

— 

— 

— 

#It ( JA-2-2-0.5 ) 

— 

— 

— 

— 

#It ( GS-1-1 ) 

80 

78 

75 

73 

#It ( GS-2-2 ) 

44 

40 

39 

44 


Table 9. Performance of the AMG solver for the linear 
elasticity problem in pure displacement form using W- 
cycle with different smoothers.. 


Level 

Li 

L 2 

L 3 


#It ( JA-1-1-0.5 ) 

50 

81 

— 

— 

#It ( JA-2-2-0.5 ) 

32 

63 

— 

— 

#It ( GS-1-1 ) 

40 

43 

43 

42 

#It ( GS-2-2 ) 

28 

27 

27 

27 


Table 10. Performance of the AMG preconditioned GG 
solver for the linear elasticity problem in pure displace¬ 
ment form using V-cycle with different smoothers. 


4.4. Numerical performance for linear elasticity problem in 
mixed form. For the linear elasticity problem in mixed displacement- 
pressure form, we plot the simulation results of the displacement and 
pressure on the left and right plots of Fig. respectively. 

We set the relative residual error l.Oe —09 in the corresponding norm 
as stopping criteria for solving the saddle point system (indehnite) with 
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Level 

Li 

L 2 


u 

#It ( JA-1-1-0.5 ) 

47 

80 

— 

— 

#It ( JA-2-2-0.5 ) 

32 

62 

— 

— 

#It ( GS-1-1 ) 

39 

44 

42 

41 

#It ( GS-2-2 ) 

28 

27 

27 

26 


Table 11. Performance of the AMG preconditioned CG 
solver for the linear elasticity problem in pure displace¬ 
ment form using W-cycle with different smoothers. 



Figure 5. Numerical results of the displacement (left) 
and pressure (right) for the linear elasticity problem in 
mixed form. 


both the AMG solver and AMG preconditioned GMRES (see [T9] l 
method. We will only consider the V-cyle for the remaining tests. 

The AMG solver with the Vanka smoother does not show the ro¬ 
bustness and efficiency in this case. However combined with GMRES 
acceleration, the V-cycle preconditioner with such a smoother shows 
improved performance. We observe acceptable performance of one or 
two V-cycles (1 V-cycle or 2 V-cycle) preconditioned GMRES solver in 
Table 12 In each cycle, we use only one pre- and post-smoothing step. 


The Braess-Sarazin smoother shows better performance. We observe 
the robustness with respect to the mesh size of the AMG solver and 
V-cycle preconditioned GMRES solver using such a smoother; see it¬ 
eration numbers of the AMG solver using one or two Braess-Sarazin 
smoothing steps (Braess-Sarazin-1-1 or Braess-Sarazin-2-2) in Table 
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Level 

Li 

L 2 

L 3 

L, 

#It ( 1 V-cycle ) 

60 

45 

49 

69 

#It ( 2 V-cycles ) 

42 

30 

32 

45 


Table 12. Performance of 


he V-cycle preconditioned 


GMRES solver for the linear elasticity problem in mixed 
form using one pre/poster Vanka smoother. 


13, and one or two V-cycles (1 V-cycle or 2 V-cycle) preconditioned 


GMRES solver in Table respectively. 


Level 

Li 

L 2 

L 3 


( Braess-Sarazin-1-1 ) 

135 

133 

125 

117 

( Braess-Sarazin-2-2 ) 

71 

70 

66 

62 


Table 13. Performance of the AMG solver for the linear 
elasticity problem in mixed form using the V-cycle with 
the Braess-Sarazin smoother. 


Level 

Li 

L 2 

L 3 


#It ( 1 V-cycle ) 

27 

27 

27 

27 

#It ( 2 V-cycles ) 

18 

19 

19 

19 


Table 14. Performance of the V-cycle preconditioned 
GMRES solver for the linear elasticity problem in mixed 
form using one pre/post Braess-Sarazin smoother. 


Using the segregated Gauss-Seidel smoother (sGS), we observe good 
performance. For all tests, we use u = 0.125. The robustness with re¬ 
spect to the mesh size of the AMG solver can be observed; see iteration 
numbers of the AMG solver using one or two segregated Gauss-Seidel 
smoothing steps (sGS-1-1 or sGS-2-2) in Table 15 The efficiency is fur¬ 
ther improved when combined with the Krylov subspace acceleration; 
see iteration numbers of one or two V-cycles (1 V-cycle or 2 V-cycle) 


preconditioned GMRES solver in Table 16 


4.5. Numerical performance for the Stokes problem. The com¬ 
putational domain for the Stokes problem is prescribed by an inside 
of a cylinder, that has radius of 1 with center point (0,0,0) on the 
inflow boundary (where u = (1.0, 0,0)), and center point (10,0,0) on 
the outflow boundary (where {2fie{u) — pl)n = (0,0,0)). On the rest 
of the boundaries u = (0, 0, 0). For all tests, we set fi = 0.5. Four levels 
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Level 

Li 

L2 

L3 


#It ( sGS-1-1 ) 

107 

104 

99 

93 

#It ( sGS-2-2 ) 

54 

53 

53 

59 


Table 15. Performance of the AMG solver for the linear 
elasticity problem in mixed form using the V-cycle with 
the segregated Gauss-Seidel smoother. 


Level 

Li 

L 2 

L 3 


#It ( 1 V-cycle ) 

14 

18 

19 

21 

#It ( 2 V-cycles ) 

12 

15 

15 

15 


Table 16. Performance of the V-cycle preconditioned 
GMRES solver for the linear elasticity problem in 
mixed form using one pre/post segregated Gauss-Seidel 
smoother. 


(Li — L 4 ) of tetrahedral meshes are generated; see mesh information 


for each level in Table 17 The number of tetrahedron (^^Tet), nodes 


(#Nodes) and midside nodes (#Midside nodes), and the total number 
of DOF ( ifDOF ) for the saddle point system. As an illustration to 
show the coarsening strategy, in Table [T^ we show the number of linear 
and quadratic velocity DOF ( 7 )^ Linear velocity DOF and # Quadratic 
velocity DOF) , the linear pressure DOF Linear pressure DOF) on 
each coarsening level for the level L 4 . The velocity and pressure of 
the simulation results are shown on the left and right plots in Fig. 
respectively. 


Level 

Li 

L2 

L3 


#Tet 

895 

7160 

57280 

458240 

#Nodes 

351 

1922 

12307 

87109 

7 )/:Midside nodes 

1571 

10385 

74802 

566212 

#DOF 

6117 

38843 

273634 

2047072 

i.BLE 17. Fluid mesh in 

brmation: Num 

Der of tet] 


hedron {^Tet) , nodes (#Nodes) and Midside nodes 
(#Midside nodes), and the total number of DOF 
(#DOF) on four levels Li — L 4 . 


For this example, the AMG solver with the Vanka, Braess-Sarazin 
and segregated Gauss-Seidel smoothers shows poor performance, that 
is very large smoothing steps are required in order to get multigrid 
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Goarsening Levels 

0 

1 

2 

3 

4 

^(/^Linear velocity DOF 

261327 

36921 

5766 

1053 

267 

T^^Quadratic velocity DOF 

1698636 

213492 

14142 

1053 

93 

pressure DOF 

87109 

12307 

1922 

351 

89 

Table 18. Number of linear anc 

quadratic velocity, anc 



linear pressure DOF in the new coarsening strategy at 
level L 4 for the Stokes problem. 



Figure 6. Numerical results of the Stokes velocity (left) 
and pressure (right). 


convergence rate. However, this will lead to very expensive computa¬ 
tional cost. In addition, we observe unsatisfactory performance of the 
AMG preconditioned Krylov subspace method using the V-cycle with 
the Vanka and segregated Gauss-Seidel smoothers. Therefore, we only 
report the performance of the AMG preconditioned GMRES solver us¬ 
ing the Braess-Sarazin smoother, that is shown in Table We set 
relative residual error l.Oe — 09 in the corresponding norm as stopping 
criteria of the AMG preconditioned GMRES solver. It is easy to see, 
with the Krylov subspace acceleration, the performance is greatly im¬ 
proved, using one or two V-cycle (1 V-cycle or 2 V-cycle) preconditioner 
with one pre- and post-smoothing steps. 
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Level 

Li 

L 2 

-^3 

L, 

#It ( 1 V-cycle ) 

25 

42 

38 

39 

#It ( 2 V-cycles ) 

16 

17 

18 

21 


Table 19. Performance of the V-cycle preconditioned 
GMRES solver for the Stokes problem using one pre/post 
Braess-Sarazin smoother. 


5. Conclusions 

In this work, we have developed an AMG method used as a stand¬ 
alone solver or preconditioner in the Krylov subspace methods for solv¬ 
ing the hnite element equations of the vector Laplacian problem, lin¬ 
ear elasticity problem in pure displacement and mixed displacement- 
pressure form, and the Stokes problem in mixed velocity-pressure form 
in 3D. We have developed a new strategy to construct the hierarchy of 
the AMG coarsening system using the hierarchical quadratic basis func¬ 
tions. The numerical studies have demonstrated the good performance 
of the AMG solvers or the AMG preconditioned Krylov subspace meth¬ 
ods for the elliptic and saddle point systems, respectively. In particular, 
the AMG preconditioned Krylov subspace methods show much better 
robustness and efficiency for solving both systems compared with the 
AMG stand-alone solvers. From this point of view, the AMG method 
developed in this work can be used as a robust and efficient solver or 
preconditioner for the SPD system and the saddle point system with 
compressible materials, and as a robust and efficient preconditioner 
for the saddle point system with incompressible materials. It is also 
possible to extend this AMG method for high-order hierarchical hnite 
element basis functions. 
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